#ifndef _h_geo
#define _h_geo

#include "bmp.h"
#include <cmath>
#include <cassert>

const double EPS=1E-12;
const double PI=acos(-1.);
const double INFTY=1E50;

int sgn(double x) {
	return x<-EPS ? -1 : x>EPS;
}

double det3(double a1, double b1, double c1,
		double a2, double b2, double c2,
		double a3, double b3, double c3) {
	return a1*b2*c3 + b1*c2*a3 + c1*a2*b3 - a1*c2*b3 - b1*a2*c3 - c1*b2*a3;
}

bool solve_e3(double a1, double b1, double c1, double d1, 
		double a2, double b2, double c2, double d2, 
		double a3, double b3, double c3, double d3,
		double *x, double *y, double *z) {
	double D=det3(a1, b1, c1, a2, b2, c2, a3, b3, c3);
	if(D==0)
		return false;
	*x=det3(d1, b1, c1, d2, b2, c2, d3, b3, c3) / D;
	*y=det3(a1, d1, c1, a2, d2, c2, a3, d3, c3) / D;
	*z=det3(a1, b1, d1, a2, b2, d2, a3, b3, d3) / D;
	return true;
}

struct Point {
	double x, y, z;
	Point()
		: x(0), y(0), z(0) {}
	Point(double x, double y, double z)
		: x(x), y(y), z(z) {}
	double get_dim(int i) const {
		assert(i>=0 && i<3);
		if(i==0) return x;
		if(i==1) return y;
		if(i==2) return z;
		return 0.;
	}
	void set_dim(int i, double v) {
		assert(i>=0 && i<3);
		if(i==0) x=v;
		if(i==1) y=v;
		if(i==2) z=v;
	}
	void uptmin_dim(int i, double v) {
		assert(i>=0 && i<3);
		if(i==0) x=min(x, v);
		if(i==1) y=min(y, v);
		if(i==2) z=min(z, v);
	}
	void uptmax_dim(int i, double v) {
		assert(i>=0 && i<3);
		if(i==0) x=max(x, v);
		if(i==1) y=max(y, v);
		if(i==2) z=max(z, v);
	}
	Point operator -() const {
		return Point(-x, -y, -z);
	}
	bool operator ==(const Point &p) const {
		return sgn(x-p.x)==0 && sgn(y-p.y)==0 && sgn(z-p.z)==0;
	}
	Point operator +(const Point &p) const {
		return Point(x+p.x, y+p.y, z+p.z);
	}
	void operator +=(const Point &p) {
		x+=p.x;
		y+=p.y;
		z+=p.z;
	}
	Point operator -(const Point &p) const {
		return Point(x-p.x, y-p.y, z-p.z);
	}
	Point operator *(double b) const {
		return Point(x*b, y*b, z*b);
	}
	Point operator /(double b) const {
		return Point(x/b, y/b, z/b);
	}
	friend Point mult(Point p1, Point p2) {
		Point p;
		p.x=p1.y*p2.z-p1.z*p2.y;
		p.y=p1.z*p2.x-p1.x*p2.z;
		p.z=p1.x*p2.y-p1.y*p2.x;
		return p;
	}
	friend Point mult(Point p0, Point p1, Point p2) {
		return mult(p1-p0, p2-p0);
	}
	friend double dot(Point p1, Point p2) {
		return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z;
	}
	friend double dot(Point p0, Point p1, Point p2) {
		return dot(p1-p0, p2-p0);
	}
	friend double dist(Point p1, Point p2) {
		return (p1-p2).len();
	}
	friend bool sameline(Point p0, Point p1, Point p2) {
		return mult(p0, p1, p2)==Point(0,0,0);
	}
	Point neg_z() const {
		return Point(x, y, -z);
	}
	double len() const {
		return sqrt(x*x+y*y+z*z);
	}
	void normalize() {
		double l=len();
		x/=l, y/=l, z/=l;
	}
	Point normalized() const {
		Point p=*this;
		p.normalize();
		return p;
	}
};

struct Ray {
	Point s, dir;
	Ray() { }
	Ray(Point s, Point dir)
		: s(s), dir(dir) { }
};

struct Patch {
	Point p[3];
	Point normal[3];
	Color c;
	Patch() { }
	Patch(Point p1, Point p2, Point p3, Point n1, Point n2, Point n3, Color c) {
		p[0]=p1;
		p[1]=p2;
		p[2]=p3;
		normal[0]=n1;
		normal[1]=n2;
		normal[2]=n3;
		this->c=c;
	}
	Patch(Point p1, Point p2, Point p3, Color c) {
		p[0]=p1;
		p[1]=p2;
		p[2]=p3;
		normal[0]=normal[1]=normal[2]=mult(p[0], p[1], p[2]).normalized();
		this->c=c;
	}
	Point get_normal(double alpha, double beta) const {
		return normal[0]*alpha + normal[1]*beta + normal[2]*(1.-alpha-beta);
	}
	Point get_point(double alpha, double beta) const {
		return p[0]*alpha + p[1]*beta + p[2]*(1.-alpha-beta);
	}
	pair<double, double> bound(int dim) const {
		pair<double, double> b=make_pair(p[0].get_dim(dim), p[0].get_dim(dim));
		for(int i=1; i<3; ++i) {
			b.first=min(b.first, p[i].get_dim(dim));
			b.second=max(b.second, p[i].get_dim(dim));
		}
		return b;
	}
	bool get_cross(Ray ray, double *ct, double *ret_alpha, double *ret_beta) const {
		double a1=p[0].x-p[2].x, b1=p[1].x-p[2].x, c1=-ray.dir.x, d1=ray.s.x-p[2].x;
		double a2=p[0].y-p[2].y, b2=p[1].y-p[2].y, c2=-ray.dir.y, d2=ray.s.y-p[2].y;
		double a3=p[0].z-p[2].z, b3=p[1].z-p[2].z, c3=-ray.dir.z, d3=ray.s.z-p[2].z;
		double alpha, beta, t;
		if(!solve_e3(a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3, &alpha, &beta, &t))
			return false;
		if(sgn(alpha)<0 || sgn(beta)<0 || sgn(alpha+beta-1.)>0)
			return false;
		if(sgn(t)<=0)
			return false;
		*ct=t;
		if(ret_alpha!=NULL) *ret_alpha=alpha;
		if(ret_beta!=NULL) *ret_beta=beta;
		return true;
	}
	bool check_cross(Ray ray) const {
		double t;
		if(!get_cross(ray, &t, NULL, NULL))
			return false;
		return sgn(t-1.)<=0;
	}
};

struct Box {
	Point s, t;
	Box() {
		s=Point(INFTY, INFTY, INFTY);
		t=-s;
	}
	void add(Point p) {
		for(int i=0; i<3; ++i) {
			s.uptmin_dim(i, p.get_dim(i));
			t.uptmax_dim(i, p.get_dim(i));
		}
	}
	double get_surface_area() {
		double a=t.x-s.x;
		double b=t.y-s.y;
		double c=t.z-s.z;
		return 2*(a*b + a*c + b*c);
	}
	bool split(int dim, double posi, Box &b1, Box &b2) {
		if(posi<s.get_dim(dim)) return false;
		if(posi>t.get_dim(dim)) return false;
		b1=*this; b1.t.set_dim(dim, posi);
		b2=*this; b2.s.set_dim(dim, posi);
		return true;
	}
	bool inside(Point p) const {
		#define CK(w) (p.w+EPS>s.w && p.w<t.w+EPS)
		return CK(x) && CK(y) && CK(z);
		#undef CK
	}
	bool check_ray(Ray ray) const {
		/*
		Point d1=s-ray.s;
		Point d2=t-ray.s;
		Point d=ray.dir;
		Point c=Point(1, 1, 1);
		if(sgn(d.x)!=0) c.y*=d.x, c.z*=d.x;
		if(sgn(d.y)!=0) c.x*=d.y, c.z*=d.y;
		if(sgn(d.z)!=0) c.x*=d.z, c.y*=d.z;
		double t1=-INFTY, t2=INFTY, p1, p2, tmp;
		for(int i=0; i<3; ++i) {
			double dx=d.get_dim(i);
			if(sgn(dx)!=0) {
				double cx=c.get_dim(i);
				p1=d1.get_dim(i)*cx;
				p2=d2.get_dim(i)*cx;
				if(p1>p2) tmp=p1, p1=p2, p2=tmp;
				if(p1>t1) t1=p1;
				if(p2<t2) t2=p2;
				if(t1>t2+EPS) return false;
			} else if(d1.get_dim(i)<-EPS || d2.get_dim(i)>EPS)
				return false;
		}
		/*/
		//*
		double t1=0, t2=1.;
		for(int i=0; i<3; ++i) {
			if(sgn(ray.dir.get_dim(i))!=0) {
				double p1=(s.get_dim(i)-ray.s.get_dim(i))/ray.dir.get_dim(i);
				double p2=(t.get_dim(i)-ray.s.get_dim(i))/ray.dir.get_dim(i);
				if(p1>p2) swap(p1, p2);
				t1=max(t1, p1);
				t2=min(t2, p2);
			} else {
				double tmp=ray.s.get_dim(i);
				if(sgn(tmp-s.get_dim(i))<0 || sgn(tmp-t.get_dim(i))>0)
					return false;
			}
			if(sgn(t1-t2)>0)
				return false;
		}
		//*/
		/*
		Point d1=s-ray.s;
		Point d2=t-ray.s;
		Point d=ray.dir;
		double t1=-INFTY, t2=INFTY;
		#define PRO(w) \
			if(sgn(d.w)!=0) { \
				double p1=d1.w/d.w; \
				double p2=d2.w/d.w; \
				if(p1>p2) swap(p1, p2); \
				if(p1>t1) t1=p1; \
				if(p2<t2) t2=p2; \
				if(t1>t2+EPS) return false; \
			} else if(d1.w<-EPS || d2.w>EPS) return false;
		PRO(x)
		PRO(y)
		PRO(z)
		//*/
			/*
		for(int i=0; i<3; ++i) {
			if(sgn(ray.dir.get_dim(i))!=0) {
				double p1=(s.get_dim(i)-ray.s.get_dim(i))/ray.dir.get_dim(i);
				double p2=(t.get_dim(i)-ray.s.get_dim(i))/ray.dir.get_dim(i);
				if(p1>p2) swap(p1, p2);
				t1=max(t1, p1);
				t2=min(t2, p2);
			} else {
				double tmp=ray.s.get_dim(i);
				if(sgn(tmp-s.get_dim(i))<0 || sgn(tmp-t.get_dim(i))>0)
					return false;
			}
			if(sgn(t1-t2)>0)
				return false;
		}
		//*/
		return true;
	}
};

#endif
